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Abstract. The subject of this paper is beam deconvolution in small angular scale CMB 
experiments. The beam effect is reversed using the Jacobi iterative method, which was de- 
\Q signed to solved systems of algebraic linear equations. The beam is a non circular one which 

£f>^ , moves according to the observational strategy. A certain realistic level of Gaussian instru 

mental noise is assumed. The method applies to small scale CMB experiments in general 
C*") . (cases A and B), but we have put particular attention on Planck mission at 100 GHz (cases 

o 

C and D). In cases B and D, where noise is present, deconvolution allows to correct the 
main beam distortion effect and recover the initial angular power spectrum up to the end of 
the fifth acoustic peak. An encouraging result whose importance is analyzed in detail. More 

5h ' 

^_ > ■ work about deconvolution in the presence of other systematics is in progress. 

<Z3 . 

■ This paper is related to the Planck LFI activities. 

> : 

• i-H , 
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1. Introduction 

Many experiments have been designed to measure Cosmic Microwave Background (CMB) 
anisotropies at small angular scales. Recent and new generation of experiments make use of 
multi-frequency and multi-beam instruments at a focus of a meter class telescope. Since not all 
the feeds can be placed along the optical axis of the telescope, the majority of them are necessary 
off-axis and the corresponding beams show more or less relevant optical aberrations (see e.g. 
Page et al. 2003 for a recent discussion of the main beam shape in the context of the WMAP 
project), according to the experiment optical design. For example, in the contex of the ESA 
Planck project 1 (see e.g. Bersanelli et al. 1996, Tauber 2000), significant efforts have been car- 
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ried out to significantly reduce the main beam distortions produced by optical aberrations (see 
e.g. Villa et al. 1998, Mandolesi et al. 2000a). On the other hand, even optimizing the optical 
design, a certain level of beam asymmetry can not be completely eliminated (see e.g. Sandri et 
al. 2002, 2003). 

If not taken into account in the data analysis, the main beam distortion introduces a system- 
atic effect in the data (Burigana et al. 1998, Mandolesi et al. 2000b) that affects the reconstructed 
map quality and, in particular, the recovery of the angular power spectrum of the CMB anisotropy 
(Burigana et al. 2001, Arnauet al. 2002, Fosalbaet al. 2002). The main topic of this paper is beam 
deconvolution in this type of experiments with the aim of remove the main beam distortion ef- 
fect in the recovery of the angular power spectrum of CMB temperature anisotropy. We wish 
to reverse the weighted average performed by a non-circular rotating beam in the presence of a 
significant level of uncorrected instrumental noise. The true angular power spectrum (Q quanti- 
ties before beam smoothing) should be recovered from the deconvolved maps, at least, in a large 
enough interval (t min , £ max ). 

A preliminary work about beam deconvolution was presented in Arnau & Saez (2000). In that 
paper, two methods for beam deconvolution were considered. One of them (hereafter method I) 
is based on the Fourier transform, and the other one (method II) uses the Jacobi algorithm for 
solving algebraic systems of linear equations. Applications of these methods in very simple cases 
were presented. The first method was applied in the case of elliptical non-rotating beams in the 
presence of a very low level of Gaussian instrumental noise. The other method was only used for 
a spherical beam in the total absence of noise. More realistic situations must be considered. This 
is the main goal of this paper. 

The formalism of our approach to deconvolution is presented in Sect. 2. 

Beam deconvolution can be only studied using simulations. Map making requires a pixeli- 
sation, and the accuracy of the angular power spectrum obtained from pixelised maps strongly 
depends on experiment sensitivity and resolution and on the sky coverage. In the case of small 
angular scales (i i 100), the angular power spectrum can be well estimated using small squared 
maps with edges lesser than 20° (Saez et al. 1996). In this case, a good map making algorithm 
and an appropriate power spectrum estimator were described in (Arnau et al. 2002). In a first 
step (first part of Sect. 3), we work with this type of simulations by assuming a simple elliptical 
main beam shape. An observational strategy involving repeated measures in the same pixel but 
without a detailed reference to a specific experiment is adopted at this stage. 

Afterwards, we apply the method to more complicated simulations carried out by using the 
HEALPix 2 (Hierarchical Equal Area and IsoLatitude Pixelization of the Sphere) package by 
Gorski et al. (1999) to pixelise the maps and compute the angular power spectrum from them. The 
beam size, its asymmetry, the variance of the instrumental noise, and the beam rotation associated 
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to an observational strategy simulating the Planck observations at 100 GHz are considered in the 
second part of Sect. 3. 

Finally, we discuss our results and draw the main conclusions in Sect. 4. 

We work in the framework of an inflationary flat universe (adiabatic fluctuations) with dark 
energy (cosmological constant) and dark matter (cold). In this ACDM model, the density pa- 
rameters corresponding to baryons, dark matter, and vacuum, are Clb = 0.05, Q.d = 0.25 and 
D.A = 0.7, respectively, and the reduced Hubble constant is h — 0.65. No reionization is consid- 
ered at all. All the simulations are based on the CMB angular power spectrum corresponding to 
this model, which has been computed with the CMBFAST code by Seljak & Zaldarriaga (1996). 

2. Beam smoothing and deconvolution 

Let us begin with an asymmetric non-rotating beam which smoothes a map T to give another 
map T s , In the continuous formalism, we can write: 



where the beam is described by function B. If the beam centre points towards a point with spher- 
ical coordinates (0, <p), the weight associated to another direction (0', <p') is a function of the form 
B{0-0',(f>-4>'). 

In the absence of rotation, function B only depends on the differences 0-0' and 0-0' 
and, consequently, Eq. ([0 is a mathematical convolution. In such a case, the Fourier transform 
can be used (as it was explained by Arnau and Saez 2000) to perform beam deconvolution. 
Nevertheless, if the asymmetrical beam rotates (as a result of the observational strategy), the 
function describing the beam is of the form B = B(0- 0' , <p — 0', 0, <f>). Since the beam is different 
(distinct orientations) when its centre points towards different points in the sky, a new dependence 
on the angles and has appeared. With a B function involving this dependence, Eq. Q is not 
a mathematical convolution and the method I, which is based on the Fourier transform, does not 
work. 

In practice, non-circular beams rotate due to the observational strategy and, moreover, the 
effect of this rotation is not negligible in many cases [see Arnau et al. 2002 for an estimation and 
Burigana et al. 2001 for an application to Planck LFI (Low Frequency Instrument, Mandolesi et 
al. 1998)]. In this situation, method I cannot be used; however, as we are going to show along the 
paper, method II works. 

Let us assume a certain pixelisation and an asymmetric beam which smoothes maps of the 
CMB sky. We first consider that only one observation per pixel is performed. The beam could 
have either the same orientation for all the pixels or different orientations for distinct pixels; in 
both cases, the beam smoothes the sky temperature T to give T s according to the relation: 




(1) 



M 




(2) 
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where M is the total number of pixels in the map. Quantity B, y is the beam weight corresponding 
to pixel j when the beam centre points inside pixel i, Eq. (|2ji can be seen as a system of M linear 
algebraic equations, in which, the independent terms T' s are the observed temperatures, and the 
M unknowns are the true sky temperatures TK The Jacobi method can be tried in order to solve 
this system. The solution would be the deconvolved map with temperatures T . No problems 
with pixel dependent beam orientations (rotation) are expected a priori; nevertheless, rotating 
asymmetric beams lead to a By matrix which is more complicated than that corresponding to a 
circular beam (studied in Arnau & Saez 2000); by this reason, we are going to verify that Jacobi 
method works for any beam, first in the absence of noise and, then, when there is an admissible 
noise level. In matrix form, Eq. @ can be written as follows: T s = BT. 

If we now consider that each pixel is observed N times either with an unique beam and 
different orientations per pixel or with various non-circular rotating beams (as it occurs in projects 
as Planck where there are various beams for each frequencies), then, we can write N matrix 
equations (one for each measure) of the form = B^ a) T, where index a ranges from 1 to N. 
The average temperature assigned to pixel i is T' a = (l/N) 2„=i T s an d the above system of N 
matrix equations leads to 



hence, for a given experiment involving various measurements per pixel, the average beam @ 
might be calculated at each pixel and, then, we can try to use the Jacobi method to solve the 
system of linear equations 0. 

We see that, in the absence of noise, method II could work in the most general case, in which 
various beams move according to the most appropriate observational strategy. For a multi-beam 
experiment one should in principle simulate the effective scanning strategy and the convolution 
with the sky signal for each beam and then apply the formalism described above by taking into 
account the various resolutions and shapes of the beams. On the other hand, the power of this 
method is that it works independently of the small differences between the resolutions and shapes 
of the various beams at the same frequency in a given experiment. Therefore, considering the 
data from a single average beam, instead of the data from the whole set of beams, but with the 
sensitivity per pixel obtained by considering the whole set of receivers at the given frequency in 
the case in which the noise is taken into account, allows to reduce the amount of data storage and 
simplify the analysis without introducing a significant loss of information about the accuracy of 
the method. 

Instrumental uncorrected Gaussian noise makes beam deconvolution more problematic; nev- 
ertheless, as we demonstrate in this work (Sects. 3.2 and 3.4), method II works for experiments 



T a = 



B a T, 



(3) 



where matrix B a describes the average beam, which can be calculated as follows: 




(4) 
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with a level of noise similar to that of Planck, or better, since the effect of deconvolution on the 
noise can be quite accurately understood with Monte Carlo simulations. 

3. Applying method II 

In this section, the deconvolution method II is applied in various cases using appropriate simula- 
tions. First, in cases A and B, method II is applied to deconvolve a set of small sky patches with 
regular pixelisation. To make this part of the work almost independent of the detailed scanning 
strategy of the considered CMB anisotropy experiment, we adopted an observational strategy 
involving multiple observations of a given pixels and only roughly mimicking that of Planck. 
The beam shape is assumed to be elliptical. Case A does not involve any noise. Case B is iden- 
tical to case A except for the presence of noise. Method II is then applied to deconvolve larger 
sky patches but using the HEALPix package for the sky pixelization and the computation of the 
angular power spectrum from coadded and deconvolved maps, simulating the Planck observa- 
tional strategy, and assuming one of the beam shapes simulated in the past year for LFI, both in 
the absence of noise (case C) and in the presence of noise (case D). 

3.1. Case A: noiseless, small patches 

An elliptical beam of the form 

f ia-»') 2 w-a') 2 1 

B(0- 6', cf> -(/>') = W N el (5) 

is assumed, as in Burigana et al. (1998) and Arnau et al. (2002). 

It is also assumed that quantities cr fl and cr^ obey the relations crgcr^ - o 2 and cr e /cr^ = 1.3, 
where cr = 4.54' (0 FWHM 10.68')- With this choice, the elliptical beam mimics the 100 GHz 
Planck beams for some locations of the detectors over the focal plane. 

We simulate squared 14.6° x 14.6° patches, with 256 (128) nodes per edge; thus, our pixel size 
is A = 3.43' (A = 6.86'). These sizes are allowed by HEALPIX and, consequently, this choice 
will facilitate some comparisons. We use seventy five of these regions covering about the forty per 
cent of the sky. With this coverage and A = 3.43' (A = 6.86'), the angular power spectrum can be 
estimated (from simulated maps) with good accuracy from €„„„ = 100to£ ma t = 10800/ A 3100 
{(max — 1550). See Saez & Arnau (1997) for details about partial coverage. In the theoretical 
model under consideration (see Sect. 1), the CMB temperature is a Gaussian homogeneous and 
isotropic statistical two dimensional field. In such a case, a certain method proposed by Bond & 
Efstathiou (1987) can be used to make the 14.6° x 14.6° maps used in this paper. This method is 
based on the following formula: 

^= 2 DVuMe-***™ , (6) 

Si,s 2 =-N 

where (\ = 2nsi/A, ti - 2ns2/A, and A stands for the angular size of the square to be mapped. 
This equation defines a Fourier transform from the position space (6, (f>) to the momentum space 
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(/i, li). The Gaussian quantities D(l\,li) have zero mean, and their variance is proportional to 
Q, where I = {t\ + (\) xl2 . Since 6T/T is real, the relation D(-t x ,-£ 2 ) = D*({ l J 2 ) must be 
satisfied. From given Q coefficients, the above D{1\ , (2) quantities can be easily calculated and, 
then, according to Eq. 0, a Fourier transform leads to the map. Saez et al. (1996) used this map 
making algorithm to get very good simulations of 20° x 20° squared regions. 

In the case of small squared maps, the above map making method suggests the power spec- 
trum estimator used in Sects. 3.1 and 3.2 and also in Arnau et al. (2002). Given one of these 
maps ST/T(0, (/>), an inverse Fourier transform leads to quantities D(l\,l2) an d, then, the average 
{\D{(i, (2)?') can be calculated on the circumference l 2 — I 2 + ^. Some interpolations are neces- 
sary to get the D(l\,ti) values at the points located on the circumference. The resulting average 
is proportional to Q , where t is the radius of the circumference. 

For A = 3.43' (A = 6.86'), the beam average can be restricted to a square with seventeen 
(nine) nodes per edge. Outside this square, beam weights given by Eq. appear to be negligible 
in this context 3 . 

Our elliptical beam is rotating while it covers a given 14.6° x 14.6° patch. In order to simulate 
beam rotation (see Fig. 1), the squared patch is located with random orientation (angle a) in the 
plane (0 , <f> ), and a different beam orientation is assigned to each pixel of the patch. If Q is 
the centre of a certain pixel, we find a point P on the 0-axis which is the centre of an auxiliary 
circumference with radius 85° which passes by Q and, then, the beam orientation - in the pixel 
under consideration - is fixed by assuming that the major axis of the elliptical beam is tangent 
- at point Q - to the auxiliary circumference. The distance from the patch centre, C, to the 8- 
axis is random, but it is constrained to be smaller than 74° in order to ensure the existence of an 
auxiliary circunference passing by the centre of every pixel. 

For each patch, the sky (T field) is simulated using either 256 or 128 nodes per edge and, 
afterwards, the beam described above is used to get the smoothed map (T s ). The pixel tempera- 
tures of the T s map are the independent terms of Eq. Q and, moreover, the terms of the B matrix 
can be built up when necessary using Eq. (jSJl and beam orientation. Taking into account that all 
the terms of the B diagonal are identical to the central weight of the beam b, the n + 1 iteration 
of the Jacobi method can be written as follows: 

r (n+l) _ T (n) + b -\j s _ b -l BT (n) _ (?) 

where T (jl) is the previous one. At zero order, we take r (0) = T s . 

Since the map is a 256x256 (128 x 128) square and the beam is another 17 x 17 (9x9) square 
(see Fig. 2), when the beam centre points towards a pixel located outside the ninth (fifth) row or 
column (counting from the nearest boundary), there are no map temperatures to be weighted. In 
practice, for CMB maps, we have verified that the following procedure works very well: Write 
3 The signal entering at angles larger than ~ 1° from the beam centre direction produces the so-called 
straylight contamination, dominated by the Galactic emission (see e.g. Burigana et al. 2001, 2003), a sys- 
tematic effect different from the main beam distortion effect considered in this work. 
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an equation for every internal node where the beam average is well defined and, then, solve 
the resulting system, which has 240 x 240 (120 x 120) equations and the same unknowns. The 
remaining temperatures (external points) are used when required by beam smoothing, but they 
are not altered along the iterative process. 

Fig. 3 shows the main results obtained in case A. Top (bottom) left panel shows quantities 
i{i+\)Cel2n in units of /iK 2 before smoothing (continuous line) and after deconvolution (pointed 
line) for A = 3.43 (A = 6.86). Both curves are indistinguishable except at the largest i values 
included in the Figure. The relative deviations between the dotted and dashed lines of each panel 
are given in the corresponding right panels. The relative error introduced by deconvolution - in 
the absence of noise - is smaller than 5% (0.5 %) for t < 1900 (t < 1480) in the case A = 3.43 
(A = 6.86); the deviations grow beyond the sixth (fifth) acoustic peak. 



3.2. Case B: noisy, small patches 

In this case, beam, rotation, coverage and pixelizations are identical to those of case A; however, 
there is instrumental uncorrelated Gaussian noise with <r N — 9 fiK (cr N = 4-.5j.iK), in antenna 
temperature, for A = 3.43 (A = 6.86), just the noise expected by combining all the beams of 
Planck at 100 GHz. A joint treatment of the impact of main beam distortions and of correlated 
1 //" type noise (see e.g. Seiffert et al. 2002) and other kind of instrumental systematics (see 
e.g. Mennella et al. 2002) is out of the scope of this paper. On the other hand, this does not 
represents a crucial limitation, since blind destriping algorithms can strongly reduce the impact 
of these effects (see e.g. Delabrouille 1998, Maino et al. 1999, Mennella et al. 2002) also in the 
presence of optical distortions (Burigana et al. 2001) and, possibly, of non negligible foreground 
fluctuations (Maino et al. 2002). The system to be solved has the form: 

M 

r s = Yj B ij Ti + N ' > ( 8 > 

!=1 

where N' is the noise at pixel i. Using matrices, this equation can be written as follows: 

T s = B(T + B~ l N) . (9) 

This last equation is formally identical to the matrix form of Eq. and it can be solved in the 
same way - using Jacobi method - to find the map T + B l N. After applying this method, some 
numerical error E is expected and, consequently, the numerical solution of system is of the 
form 

T* = T + B~ l N + E. (10) 

In general, T* is different from T (sky temperature before smoothing); hence, the angular power 
spectrum extracted from the map T* is different from that of the unconvolved sky, which is to be 
extracted from map T. Results are shown in Fig. 4, which has the same structure as Fig. 3. We see 
that the spectra before smoothing and after deconvolution (which are obtained from maps T and 
T*, respectively) separate at middle t values. In the right panels, we can verify that the deviation 
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produced by deconvolution - in the presence of the assumed level of noise - is smaller than five 
per cent for I < 1 100 (I < 950) in the case A = 3.43 (A = 6.86). The relative deviations rapidly 
increases for greater I values. From the comparison of Figs. 3 and 4 it follows that the presence 
of noise with cr N = 9/iK (<r N = 4.5/i/f) has important consequences for beam deconvolution. 

Fortunately, the angular power spectrum of the map N* = B~ l N can be estimated and sub- 
tracted from that of the T* map. In order to do that, we first solve the matrix equation: 

N' = BN*. (11) 

This equation can be also solved using the Jacobi method. The independent terms are the temper- 
atures N' corresponding to a new noise realization different from N. We can now extract the Q 
quantities from the map N*. New numerical errors should appear when we use the Jacobi method 
in the presence of noise. On the other hand, we can try a Monte Carlo approach to evaluate the 
deconvolution effect on the noise. We can take various N' noise realizations to get an average 
spectrum of the corresponding N maps; in this way, the effect of the noise variance in the esti- 
mate of the angular power spectrum of N* is strongly reduced (forty noise realizations suffice). 
When we subtract this spectrum from that of T", namely, when we correct the T* spectrum tak- 
ing into account noise effects, results are much better than those showed in Fig. 4 4 . These new 
results are presented in Fig. 5. The structure of this figure is identical to that of Figs. 3 and 4. As 
in Fig. 3, the range of I values - in the left panels - has been appropriately chosen to include the 
region where the displayed curves separate significantly. The deconvolution is better than that of 
Fig. 4, where no correction for the noise has been considered. According to the right panels of 
Fig. 5, the relative deviation produced by deconvolution plus correction is smaller than five per 
cent for I < 1500 (£ < 1300) in the case A = 3.43 (A = 6.86). For A = 3.43 (A = 6.86), decon- 
volution works very well up to the end of the fifth (fourth) acoustic peak. For equivalent levels 
of noise in different pixels, we see that - as expected - deconvolution has recovered more Q 
coefficients for A = 3.43; hence, we can say that results corresponding to A = 3.43 are sensibly 
better than those of A = 6.86. On account of this fact and also for the sake of briefness, pixels 
with A = 6.86 are not considered in cases C and D below. 

3.3. Case C: application to Planck in the absence of noise 

The selected orbit for Planck is a Lissajous orbit around the Lagrangian point L2 of the Sun- 
Earth system (see e.g. Bersanelli et al. 1996). The spacecraft spins at 1 r.p.m. and the field of 
view of the two instruments - LFI and HFI (High Frequency Instrument, Puget et al 1998) - 
is about 10° X 10° centered at the telescope optical axis (the so-called telescope line of sight, 
LOS) at a given angle a from the spin-axis direction, given by a unit vector, s, chosen to be 
4 We observe that an analogous approach can be pursued also in the presence of correlated noise, provided 
that the noise properties can be known from laboratory measures and/or directly reconstructed from the data 
(Natoli et al. 2002). Of course, in this context, destriping (or, possibly, map-making, see e.g. Natoli et al. 
2001) should be previously applied both to the data and to the simulated pure noise data. 
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pointed in the opposite direction with respect to the Sun. In this work we consider values of 
a ~ 85°, as adopted for the baseline scanning strategy. The spin axis will be kept parallel to the 
Sun-spacecraft direction and repointed by 2.5' every 1 hour (baseline scanning strategy). 
Hence Planck will trace large circles in the sky and we assume here, for simplicity, 60 exact 
repetitions of the set of the pointing directions of each scan circle. A precession of the spin-axis 
with a period, P, of 6 months at a given angle ft ~ 10° about an axis, /, parallel to the Sun- 
spacecraft direction (and outward the Sun) and shifted of 2.5' every 1 hour, may be included 
in the scanning strategy, possibly with a modulation of the speed of the precession in order to 
optimize data transmission (Bernard et al. 2002). The quality of our deconvolution code is of 
course almost independent of the details of these proposed scanning strategies, and we assume 
here the baseline scanning strategy for sake of simplicity. 

The code implemented for simulating Planck observations for a wide set of scanning strate- 
gies is described in detail in Burigana et al. (1997, 1998) and in Maino et al. (1999). In this study 
we do not include the effects introduced by the Planck orbit, to be currently optimized, by simply 
assuming Planck located in L2, because they are fully negligible in this context. 

We compute the convolutions between the antenna pattern response and the sky signal as 
described in Burigana et al. (2001) by working at ~ 3.43' resolution and by considering spin- 
axis shifts of ~ 2.5' every hour and 7200 samplings per scan circle. We simulate 11000 hours 
of observations (about 15 months) necessary to complete two sky surveys with the all Planck 
beams. 

With respect to the reference frames described in Burigana et al. (2001), following the recent 
developments in optimizing the polarization properties of LFI main beams (see e.g. Sandri et al. 
2003), the conversion between the standard Cartesian telescope frame x, y, z and the beam frame 
Xbf,ybf,Zbf actually requires a further angle tfrs other than the standard polar coordinates 8b and 
(f>B defining the colatitude and the longitude of the main beam centre direction in the telescope 
frame. Appendix A provides the transformation rules between the telescope frame and the beam 
frame, as well as the definition of the reference frames adopted in this work. 

The orientation of these frames as the satellite moves is implemented in the code. For each 
integration time, we determine the orientations in the sky of the telescope frame and of the beam 
frame, thus performing a direct convolution with the sky signal by exploiting the detailed main 
beam response in each considered sky direction. The detailed main beam shape and position on 
the telescope field of view adopted in this application is that computed in the past year for the feed 
LFI9 (Sandri et al. 2002) which shows an effective FWHM resolution of 10.68' and deviations 
from the symmetry producing a typical ellipticity ratio of 1.25. Such values of resolution and 
asymmetry parameter are in the range of them that it is possible to reach with a 1 .5 m telescope 
like that of Planck by optimizing the optical design (see e.g. Sandri et al. 2003). Although our 
deconvolution method is largely independent of the details of the considered beam shape, it is 
interesting to exploit its reliability under quite realistic conditions. 



10 C. Burigana & D. Saez: Beam deconvolution in noisy CMB maps 

The CMB anisotropy map has been projected into the HEALPix scheme (Gorski et al. 1999) 
starting from the angular power spectrum of the assumed ACDM model (see Sect. 1). 

To make the application of the deconvolution code easier and the system solution possible 
without large RAM requirement and in a reasonable computational time 5 we implemented a 
code that identify in the simulated time ordered data (TOD) all and only the beam centre pointing 
directions in an equatorial patch (in ecliptic coordinates) of 1024 x 1024 pixels with a ~ 3.43' 
side (n S ide = 1024). We keep the exact information on the beam centre pointing direction and 
the beam orientation (defined for instance by an angle between the axis x\,f and the parallel in 
the beam centre pointing direction) as computed by our flight simulator. All the samples of the 
TOD within the same pixel are identified and restored in contiguous positions. At this aim, we 
take advantage from the nested, hierarchical ordering of the HEALPix. This is quite simple in the 
current simplified simulation, but it will require the development of efficient and versatile tools 
to manage the more general case in which the all samples from the experiment multi-beam array 
are considered, particularly for the ecliptic polar patches, which pixels are observed many and 
many times because of the Planck scanning strategy. In the context of the Planck project, this 
effort will be pursued by taking advantage from the development of Planck Data Model (see e.g. 
Lama et al. 2003). 

^From the simulated TOD, possibly restored as described above, we extract a map of a patch 
of simply coadded data and a map of a patch deconvolved by applying method II. The latter 
map can be then symmetrically smoothed with a beam FWHM of 10.68' by using the HEALPix 
tools for comparison with the former one, obtained from the convolution with the simulated 
asymmetric beam and taking into account the scanning strategy. Of course, from the input map 
we can extract the same sky patch. 

We consider four different patches covering an equatorial region of ^ 28.3 % of the sky 
(analogously to the case of small patches, see Sect. 3.1, avoiding the boundary regions of the 
four patches slightly reduces the originally considered, ^ 33.3 %, sky coverage). 

All the above maps are inverted with the anafast code of HEALPix to extract the corre- 
ponding angular power spectra. The result is shown in Fig. 6. Of course, all the angular power 
spectra are in strict agreement at multipoles t 200, where the main beam distortion effect is 
negligible for a beam with a FWHM of about 10' and a reasonable ellipticity. Note the very 
good agreement between the power spectrum of the input map and that derived with method II: 
the difference becomes significant only at the seventh acoustic peak [compare the solid (black) 
line with the dashed (green) line]. Note the power excess at high ( introduced by the beam dis- 
tortions, compared with the power spectrum derived from the deconvolved map subsequently 
symmetrically smoothed [compare the dash-three dots (fuxia) with the dash-dots (blue)]. Also, 
the power spectrum derived from coadded map when divided by the window function corre- 

5 In the current implementation, about 19 hours of computation are required to deconvolve a single patch 
with 1024 2 pixels on an 64 bit alpha digital unix machine with single cpu at 533 MHz and 1 Gb RAM. 
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sponding to the symmetric equivalent beam, exp[-(cr{) 2 ], significantly exceeds that of the input 
map at t larger than the fourth acoustic peak [compare the dots (red) with the solid line (black)]. 
We find a similar disagreement even by varying the assumed value of the symmetric beam width 
<x: an improvement on a limited range of multipoles results in a worsening on a different range 
of multipoles. 

This demonstrates that a kind of deconvolution is necessary to remove the main beam dis- 
tortion effect at very high multipoles and that method II works well in the absence of noise. The 
impact of instrumental noise is discussed in the next subsection. 

3.4. Case D: application to Planck in the presence of noise 

Analogously to the Case B (see Sect. 3.2) we have simulated four patches (of 1024 2 pixels, as 
in the previous section) of instrumental uncorrected Gaussian noise with ctn = 9//K (in terms 
of antenna temperature) for a pixel of 3.43', as appropriate to the global Planck sensitivity at 
100 GHz; for simplicity we have assumed a uniform noise. 

This realization of noise map has been added to the map of signal and then method II has been 
applied to deconvolve the map of signal plus noise, as described in Sect. 3.2. We produced also 
four realizations of pure noise maps to be deconvolved in the same way. Finally, we generated 
another map of noise to be superimposed to the coadded map obtained from the convolution with 
the simulated beam including scanning strategy, for comparison. 

We computed the angular power spectrum of the four maps of pure noise and of the 
four maps of pure noise deconvolved with method II. Fig. 7 (left panel) compares the aver- 
ages of the four realizations of these power spectra and their relative variance (right panel). 
As evident, deconvolution increases the noise: a rough approximation of the ratio between of 
the noise angular power spectrum after deconvolution and before deconvolution is given by 
~ 2(FWHM/A) 2 exp[(c^/2) 2 + (cr^/2) 6 ] where, as usual, FWHM = V8Tn2<x (= 10.68') and 
A is the pixel side (= 3.43'). On the other hand, the relative variance of these power spectra is 
almost similar. 

In spite of the relatively large increase of the noise power, we find that method II results to 
work quite well in removing the effect of main beam distortions up to the end of the fifth acoustic 
peak, when the average power spectrum of the deconvolved pure noise maps is subtracted to the 
power spectrum of the deconvolved noisy map (see the dashed (green) line in Fig. 8). 

In Fig. 9 we report the relative (per cent) errors introduced by beam distortions in the ab- 
sence of deconvolution, in the presence of deconvolution without applying the subtraction of 
the average deconvolved noise spectrum and by applying the deconvolution and the subtraction 
of average deconvolved noise spectrum. As evident, in the last case the power spectrum can be 
recovered with a good accuracy up to high multipoles (relative errors £ 5, 10, 15, 20 % respec- 
tively for £ i 1250, 1470, 1500, 1650 - see the dashed (green) line in the middle panel - to be 
compared with errors ^ 10, 20, 30 % at ( ^ 1200, 1250, 1300 - see the dotted (red) line in the 
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middle panel - and then dramatically increasing with t in the absence of deconvolution). The 
right panel shows what already found for the noiseless case (Sect. 3.3): even by varying the as- 
sumed value of the symmetric beam width cr, an improvement in the Q recovery can not be 
reached simultaneously on the whole relevant range of multipoles. 

We have considered here for simplicity only equatorial patches. On the other hand, we have 
verified that the main beam distortion affects the reconstructed power spectrum in a similar way 
also on polar patches. In fact, even at high ecliptic latitudes a given pixel is preferentially ob- 
served with a limited number of beam orientations. In polar patches the Planck sensitivity is 
significantly better than the average. Therefore, in spite of the more complex data storage, de- 
convolution will be there less affected by the noise 6 . 

4. Discussion and conclusions 

We have presented the basic formalism for a robust and feasible method for deconvolution in 
noisy CMB maps. 

We have implemented and tested this method for two completely different situations: small 
patches observed with an elliptical beam and with a scanning strategy involving repeated mea- 
sures of the same pixel but not related to a specific experiment and quite large sky areas observed 
with a realistic beam and a scanning strategy like that of the Planck satellite. A sensitivity level 
and a beam resolution typical of the Planck experiment at 100 GHz have been exploited. 

After having verified the good accuracy of our deconvolution code to remove the main beam 
distortion effect in the angular power spectrum estimation in the absence of noise, we applied it to 
noisy maps. We demonstrate that it is possible to accurately evaluate the effect of deconvolution 
on pure noise simulated maps so deriving, with Monte Carlo simulations, a good estimation of 
the average deconvolved noise angular power spectrum to be subtracted from the deconvolved 
noisy maps. 

In practice, for the considered sensitivity, 9yuK for a pixel of 3.43', and beam resolution, 
FWHM =! 10'- 11 ', our deconvolution code allows to efficiently remove the main beam distortion 
effect and accurately reconstruct the CMB power spectrum up to the end of the fifth acoustic 
peak. 

Clearly, in the context of the Planck project, the measure of the very high multipole region of 
the CMB angular power spectrum will be take advantage from the cosmological frequency chan- 
nels at highest resolution, namely the 217 GHz channel (having a FWHM 5'), where Poisson 
fluctuations from extragalactic sources (see e.g. Toffolatti et al. 1998) are expected to be at a 
very low level and anisotropies from thermal Sunyaev-Zeldovich effects are, if not exactly null 
because of possible unbalanced contributions within the bandwidth, certainly very small. On the 

6 For example, the average sensitivity on a polar region of about 25 squared degrees is about 5 times better 
(i.e. =! 2//K on pixel of 3.43'!) than the average full sky sensitivity: then, we expect there a deconvolution 
quality intermediate between that found here and that found in the previous section. 
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other hand, the frequency range about 100 GHz is where the global (Galactic plus extragalac- 
tic) foreground contamination is expected to be minimum. Therefore, it is extremely relevant to 
extract at these frequency an accurate estimation of the sky angular power spectrum, cleaned, 
as better as possible, from the all systematic effects. In addition, the removal of the main beam 
distortion effect, relevant at large multipoles, greatly helps the comparison between the results 
obtained at different frequency channels. 

Of course, we plan to apply this method in the future also to lowest and highest beam reso- 
lutions and in the presence of other kinds of systematic effects. We believe that the results found 
here are very encouraging, suggesting that the main beam distortion effect, previously reduced 
by optimizing the optical design, can be further reduced in the data analysis. 
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Appendix A: Transformation rules between Planck telescope frame and beam 
frame 

Let s be the unit vector, choosen outward the Sun direction, of the spin axis direction and k that 
of the direction, z, of the telescope line of sight (LOS), pointing at an angle a ~ 85° from the 
direction of s. 

On the plane tangent to the celestial sphere in the direction of the LOS we choose two coor- 
dinates x and y, respectively defined by the unit vector i and j according to the convention that 
the unit vector i points always toward s and that x, y, z is a standard Cartesian frame, referred 
here as telescope frame. 

Let it,f,jbf,kbf be the unit vectors corresponding to the Cartesian axes Xbf,ybf,Zbf of the 
beam frame; kbf defines the direction of the beam centre axis in the telescope frame. The beam 
frame is defined with respect to the telescope frame by three angles: 8b, <Pb, i'B (9b an d <Pb, two 
standard polar coordinates defining the direction of the beam centre axis, range respectively from 
0°, for an on-axis beam, to some degrees, for LFI off-axis beams, and from 0° to 360°). 

Let i' bf, j' b f,k'bf (k'bf = kbf) be the unit vectors corresponding to the Cartesian axes x'.y'.z' 
of an intermediate frame, defined by the two angles Og and <pB, obtained by the telescope frame 
x, y, z when the unit vector of the axis z is rotated by an angle Og on the plane defined by the unit 
vector of the axis z and the unit vector kbf up to reach kbf 

k'bf = kbf = cos(0 B )sin(0 B ); + sin(0 B )sin(# B ) j + cos(6 B )k (A.l) 
f'bf = [cos(0 B ) 2 cos(# B ) + sm((f> B ) 2 ]l + [sin(0 B )cos(0 B )(cos(0 B ) - l)]j - sin(0 B )cos(0 B )£ (A.2) 
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j' bf = [sin(<^)cos(</> B )(cos(6» B ) - l)]z + [cos(# B )sin(0 B ) 2 + cos(0 B ) 2 ] j - sm(0 B )sm{4> B )k . (A3) 

The beam frame is obtained from the intermediate frame through a further (anti-clockwise) 
rotation of an angle iff B (ranging from 0° to 360° 7 ) around k b f and is therefore explicitely given 
by: 



i bf = [cos((/r B >'V >A . + sm.(i/fB)j' bfx \i + [cos(i(/ B )i' bLy + sm(\j/ B )j' bfy \j 
+[cos(iff B )i' bLz + sin(tff B )j' bfy ]z 



(A.4) 



j bf = [-sm(i// B )i' bfyX + cos(i// B )j' bf Ji + [-sin(iff B )i' b fy + cos(^ fi )/' i/v ] j 

+ [sin(tl/ B )i' bf, z + cos(^b) j' bftZ ]Z , (A.5) 

where the bottom index x (y, z) indicates the component of intermediate frame unit vector along 
the axis x (y,z) of the telescope frame, as defined by eqs. (A1-A3). 
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Fig. 1. Patch location and beam orientation in cases A and B denned in the text. Points C and 
Q are the centres of the patch and the pixel, respectively. Angle a and the angular distance 
DC < 74° define patch location. Point P is the centre of a circumference with an angular radius 
of 85°, which passes by point Q. The tangent to this curve at Q defines beam orientation (see also 
the text). 
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Fig. 2. Boundary conditions for the application of the Jacobi method in cases A and B (see text). 
An equation is written at each internal node (diamonds). No equations are associated to external 
points (asterisks). Temperatures of the external nodes keep unaltered along the iterative process. 




Fig. 3. Top left panel displays quantities ((( + l)C(/2n, in units of /iK 2 , extracted from simulated 
maps with a pixel size A = 3.43. solid (dotted) line has been obtained from unconvolved (decon- 
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Fig. 4. The same as in Fig. 3 but in the presence of a level of noise cr N — 9 fiK for the pixel size 
A = 3.43 (top panels) and cr N = 4.5 fiK for A — 6.86 (bottom panels). No correction for noise 




Fig. 5. The same as in Fig. 4, but with correction for noise effects on the angular power spectrum 
(see also the text). 
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Fig. 6. The two panels are identical except for the binning in ( in the right panel. No noise is 
considered. Angular power spectrum from the four considered patches. Solid (black) line: angu- 
lar power spectrum of the input map without convolution; dash-three dots (fuxia): angular power 
spectrum of the map convolved with the simulated beam assuming the Planck scanning strategy; 
dots (red): angular power spectrum of the map convolved with the simulated beam assuming the 
Planck scanning strategy and divided by the symmetric beam window function; dashes (green): 
angular power spectrum of the deconvolved map; dash-dots (blue): angular power spectrum of 
the deconvolved map subsequently convolved with a symmetric beam with the effective resolu- 
tion (see also the text). 
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Fig. 7. Left panel: comparison between the noise power spectrum before (red) and after (green) 
deconvolution for the four noise realizations (the inner (black) "curves" represent the two average 
noise power spectra). Right panel: ratio between the power spectrum of each of the four noise 
realizations and the average noise power spectrum, before and after deconvolution, multiplied by 
10 in the latter case for graphic purposes (see also the text). 
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Fig. 8. The two panels are identical except for the binning in £ in the right panel. The noise is 
included. Angular power spectrum from the four considered patches. Solid (black) line: angular 
power spectrum of the input map without convolution and without noise; dash-three dots (fuxia): 
angular power spectrum of the map convolved with the simulated beam assuming the Planck 
scanning strategy and adding a noise realization; dots (red): angular power spectrum of the map 
convolved with the simulated beam assuming the Planck scanning strategy and adding a noise 
realization, after the subtraction of the averaged power spectrum of four noise realizations and 
then divided by the symmetric beam window function; dotted (black) bottom line: averaged 
power spectrum of four noise realizations; dash-dots (green): angular power spectrum of the 
deconvolved map in the presence of a noise realization; dashes (green): angular power spectrum 
of the deconvolved map in the presence of a noise realization after the subtraction of the averaged 
power spectrum of four noise realizations deconvolved in the same way; dashed (black) bottom 
(at low t) line: averaged power spectrum of four noise realizations deconvolved in the same way 
(see also the text). 
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Fig. 9. The left and middle panels are identical except for the binning in t in the middle panel. 
Errors in the angular power spectrum recovery from the four considered patches. Dash-three 
dots (fuxia): angular power spectrum of the map convolved with the simulated beam assuming 
the Planck scanning strategy and adding a noise realization; dots (red): angular power spec- 
trum of the map convolved with the simulated beam assuming the Planck scanning strategy and 
adding a noise realization, after the subtraction of the averaged power spectrum of four noise 
realizations and then divided by the symmetric beam window function; dash-dots (blue): angular 
power spectrum of the deconvolved map in the presence of a noise realization; dashes (green): 
angular power spectrum of the deconvolved map in the presence of a noise realization after the 
subtraction of the averaged power spectrum of four noise realizations deconvolved in the same 
way. In the right panel we report the same line (dots, red) of the middle panel and four other 
power spectra obtained in similar way by assuming different values of effective angular resolu- 
tion, with an effective FWHM increased or decreased by 0.05' and 0.1'. Clearly, an improvement 
at i ~ 400 - 800 implies a worsening at t i 1000 and viceversa: therefore, even allowing for 
changes in the assumed effective angular resolution, the simple symmetric beam approximation 
can not improve the power spectrum recovery simultaneously in the two above ranges of ( (see 
also the text). 



